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We analyze the properties of particles trapped in three-dimensional potentials formed from su¬ 
perimposed Gaussian beams, fully taking into account effects of potential anharmonicity and non¬ 
separability. Although these effects are negligible in more conventional optical lattice experiments, 
they are essential for emerging ultracold atom developments. We focus in particular on two poten¬ 
tials utilized in current ultracold atom experiments: arrays of tightly focused optical tweezers and 
a one-dimensional optical lattice with transverse Gaussian confinement and highly excited trans¬ 
verse modes. Our main numerical tools are discrete variable representations (DVRs), which combine 
many favorable features of spectral and grid-based methods, such as the computational advantage of 
exponential convergence and the convenience of an analytical representation of Hamiltonian matrix 
elements. Optimizations, such as symmetry adaptations and variational methods built on top of 
DVR methods, are presented and their convergence properties discussed. We also present a quan¬ 
titative analysis of the degree of non-separability of eigenstates, borrowing ideas from the theory of 
matrix product states (MPSs), leading to both conceptual and computational gains. Beyond devel¬ 
oping numerical methodologies, we present results for construction of optimally localized Wannier 
functions and tunneling and interaction matrix elements in optical lattices and tweezers relevant for 
constructing effective models for many-body physics. 


I. INTRODUCTION 

Ultracold neutral atoms trapped in optical potentials 
have been solidly established as a highly controllable 
platform for precision measurement and quantum metrol¬ 
ogy m as well as quantum simulation of many-body 
physics mm- The most prevalent method for quantum 
simulation of three-dimensional (3D) lattice systems is 
to optically trap neutral atoms in a cubic lattice formed 
by three orthogonal pairs of counterpropagating laser 
beams [4]. An emerging alternative for manipulating 
and trapping ultracold atoms is optical tweezers 0 , in 
which a tightly focused Gaussian beam is used to trap 
single atoms. In contrast to optical lattice experiments, 
which typically load gases which have been evaporatively 
cooled, atoms in optical tweezers can be individually ma¬ 
nipulated and cooled to their ground state using laser 
cooling alone [H]. The ability to dynamically vary the 
position of such tweezers at separations where tunnel¬ 
ing is appreciable, as has been demonstrated in recent 
experiments mm, leads to a “bottom-up” approach to 
building low-entropy quantum systems, in contrast to the 
“top-down” approach of most optical lattice experiments 
in which large ensembles of atoms are cooled in an exter¬ 
nal trap before being loaded into the lattice. 

The many-body physics of ultracold atoms in optical 
potentials is usually described within the framework of 
Hubbard models, which are truncated expansions of the 
full Hamiltonian in a basis of spatially localized orbitals 
known as Wannier functions [S]. The parameters ap¬ 
pearing in such models, for example tunneling integrals 
and interaction matrix elements, are the point of con¬ 
nection between few-body physics in the confining po¬ 
tential and many-body physics: these parameters con¬ 


trolling the many-body physics can be determined from 
few-body calculations or experiments. Hence, determin¬ 
ing them quantitatively is key to designing and vali¬ 
dating robust quantum simulators and other quantum 
technologies. The construction of Hubbard models for 
cubic optical lattices is greatly facilitated by the fact 
that such lattices are separable in Cartesian coordinates, 
U(r) = V{x) + V{y) + V{z), and hence their eigenfunc¬ 
tions are products of eigenfunctions of one-dimensional 
(ID) problems, which are much less computationally de¬ 
manding than 3D problems. In addition, the theory of 
periodic potentials in ID is especially simple and leads 
to computationally efficient procedures for determining 
the Wannier functions with maximum spatial localiza¬ 
tion m- On the other hand, technologies which depend 
upon the curvature of light beams for trapping, such 
as optical tweezers, are inherently non-separable, and so 
the well-developed machinery employed for cubic lattices 
is not applicable. In addition to being more computa¬ 
tionally challenging, non-separable potentials have sig¬ 
nificant qualitative differences from separable potentials. 
For example, the tunneling rate of a particle through a 
non-separable lattice depends on its transverse motional 
state, and can change this rate by an order of magnitude 
or more. In contrast, for a separable potential tunneling 
rates are independent of the transverse motional state. 

In this work, we study the eigenstates of non-separable 
3D optical potentials constructed by superimposing 
Gaussian laser beams, taking as our two main examples 
arrays of optical tweezers [3 [ 5 ] and a ID optical lattice 
with transverse Gaussian confinement, such as is utilized 
for optical lattice clocks m- Our main numerical tool 
is discrete variable representations (DVRs), coupled with 
variational methods. As will be discussed in Sec. 
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given as appendices. 


II. OPTICAL POTENTIALS FOR NEUTRAL 
ATOMS 

In this work, we consider optical potentials which are 
generated by superimposed Gaussian laser beams. The 
electric held amplitude of the fundamental TEMqo mode 
of a Gaussian laser beam propagating along z may be 
written as [13] 


FIG. 1: Double-well optical tweezer. The double-well optical 
tweezer potential, Eq. a,t y — 0, with no bias (A = 0) and 
spacing a/wo = 1.23. Note the different scaling of the x and 
z coordinates. 
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DVR methods combine many nice features of grid-based 
methods, such as analytic representation of the kinetic 
energy operator and a diagonal representation of the po¬ 
tential energy operator, with features of spectral meth¬ 
ods, such as exponential convergence. We focus in par¬ 
ticular on the connections between the trap parameters, 
e.g., the trap depth, and the parameters appearing in ef¬ 
fective many-body models, such as tunneling and inter¬ 
action matrix elements. We also quantitatively study the 
degree of separability of the eigenfunctions of such poten¬ 
tials using tools from quantum information science. In 
particular, borrowing from the theory of matrix product 
states (MPSs) [l2], which accurately capture many-body 
states with restricted entanglement, we develop an analo¬ 
gous state ansatz for non-separable single-particle states. 
As we will show, this ansatz is useful for visualization, 
storage, and computational efficiency. 

This work is organized as follows: Sec. [IT] defines the 
potentials we study in this work, their symmetries, and 
the Hubbard parameters we will study; Sec. |III| reviews 
the two DVR basis sets we use and presents numerical op¬ 
timizations for DVR-based algorithms; Sec. |IV| presents 
an analysis of non-separable single-particle states from 
the viewpoint of quantum information theory, in partic¬ 
ular discussing a canonical form for non-separable states 
motivated by matrix product states; Sec. [V| gives results 
for Hubbard parameters and quantihes non-separability 
of states in a double-well optical tweezer array and a non- 
separable optical lattice; Finally, Sec. |Vl] concludes and 
gives an outlook. Some technical details on computing 
interaction matrix elements in a basis of radial functions 
and a variational algorithm for finding the nearest sepa¬ 
rable state given a state in the MPS canonical form are 


where Eq is the field amplitude, the beam waist wq is 
where the field amplitude drops to 1/e of its on-axis 
value, and the Rayleigh range zn = ttWq/X, with A the 
wavelength of the laser light. The resulting optical po¬ 
tential is proportional to the intensity of the field and 
the atomic polarizability. For a single held of the form 
Fq. 0 with the laser-frequency far detuned from an opti¬ 
cal transition, the ac-Stark shift gives rise to the potential 



with maximum depth at the origin V. We assume the 
laser is red-detuned, and therefore that V is positive. 
We will refer to a single potential of the form Fq. 0 as 
an optical tweezer when the beam waist is comparable to 
the wavelength wq ~ A. Such a potential can be formed 
by focusing a Gaussian beam through a high numerical 
aperture lens. Expanding Eq. ([^ to lowest order in r 
and z results in a harmonic approximation 

Vs (r) « -V -k -k , (3) 

Wo Zr 

corresponding to radial and axial frequencies fiWr = 
yj8Eisg V and Hujz = a/TEZ/V, respectively, where 
Ewo = fr^/imwQ and = h‘^/2mzj^ are the charac¬ 
teristic energies associated with the waist and Rayleigh 
range, respectively. 

The second case we consider is that of two focused 
traps in which one has an intensity greater by A more 
than the other, resulting in a double-well optical tweezer 
potential of the form 


Vd (r) = - 



V exp — 



-k (H + A) exp - 



( 4 ) 
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FIG. 2: Non-separable optical lattice. The non-separable op¬ 
tical lattice, Eq. ([^, for y/I4onst = 3/22. A non-zero value of 
Konst increases the effective transverse (r) confinement while 
leaving the axial ( 2 ) motion unchanged, to lowest order. Note 
the different scalings of the r and 2 coordinates, as well as the 
difference in scaling of 2 coordinates with respect to Fig. [T]d ue 
to the additional length scale a. 


as shown in Fig. when the lasers are set up to avoid 
coherent interference effects, for example by frequency 
detuning the two beams with a splitting above any 
trap time scales. The dynamics of two bosonic [7] or 
fermionic [5] atoms in such a double-well tweezer config¬ 
uration have been investigated in recent experiments. 

The third case that we consider is a ID optical lattice 
with transverse Gaussian character, which can be created 
by interfering two counter-propagating beams of the form 
Eq. Q. Confinement along the 2 direction is provided 
by the standing wave interference pattern, and sowqc^X 
is not required for 3D confinement. Hence, we consider 
the regime Wq 3> A in which most experiments operate. 
In this regime, we can neglect any effects of the Rayleigh 
range, giving rise to the potential 

V (r,z) = -Kattexp cos^ (kz) , (5) 

V ■^o / 

where k = tt/o, with a = X/2 the lattice spacing. In 
addition, to study the effects of non-separable transverse 
confinement on the axial motion of particles through a 
ID optical lattice, we will consider the generalized ID 
optical lattice 



which is comprised of a standing wave optical lattice 
Eq. (§ together with the potential resulting from a non- 
reflected beam Eq. ([^ of the same waist. In this con¬ 
figuration, the transverse confinement frequency is set 
by ^/8EyJg (Konst + ^)) while confinement along 2 mea¬ 
sured close to a single lattice minimum is set by V alone. 

Atoms loaded into an optical tweezer potential can be 
efficiently cooled to their 3D ground state via laser cool¬ 
ing 0, and so our analysis of tweezer potentials will be 


focused on the properties of the lowest few eigenstates. 
In contrast, some applications of optical lattices, such as 
optical atomic clocks do not require being in the 

lowest transverse motional states. Hence, in discussing 
the ID optical lattice potential we will devote special 
attention to develop methods which can efficiently deal 
with a large number of states to facilitate thermal aver¬ 
ages. 


A. Symmetries of the optical potentials and 
characterization of eigenstates 


For a periodic potential such as the ID optical lattice 
potentials Eq. the solutions may be character¬ 

ized by a quasimomentum q in the hrst Brillouin Zone 
(BZ), q S [—tt/o, 7r/a), even in the case that the poten¬ 
tial is non-separable. The resulting eigenfunctions can be 
written in Bloch form i/nq (r) = e^'^^Unq (r) /VL, where 
Unq{x,y,z + a) = Unq{x,y,z) is a unit periodic func¬ 
tion in 2 , with a the lattice spacing and L the number 
of unit cells. Eigenfunctions of different q are orthog¬ 
onal by construction, and so can be obtained in sepa¬ 
rate calculations. The double-well potential Eq. Q has 
a mirror symmetry along the y and 2 directions. As a 
consequence, the solutions can be characterized in terms 
of their parities Py and where P„ = 1 (—1) corre¬ 
sponds to a function which is even (odd) under inversion 
of coordinate v. In addition, for A = 0 there is also 
mirror symmetry along x. As with the quasimomentum, 
we can obtain the eigenstates in each symmetry sector 
separately, leading to significant computational gains. 


For a non-separable potential, we cannot unambigu¬ 
ously assign quantum numbers counting the number of 
quanta of excitation along each direction. However, in 
most cases it is possible to assign labels which specify 
that the separable state nearest to the given state may be 
labeled by a vector of excitation quanta n = (n^, n^) for 
problems of cylindrical symmetry and n = {nx,ny,nz) 
for problems without cylindrical symmetry. This vector 
n can be determined by counting the number of nodes in 
the components of the nearest separable wavefunction. 
A precise definition of the nearest separable state and a 


procedure for obtaining it are presented later in Sec. IV 


In summary, the eigenstates of the optical lattice po¬ 
tentials Eqs. and § can be written as tpmnrn^q (r), 
where q is the quasimomentum, m is the azimuthal quan¬ 
tum number arising from cylindrical symmetry, and n^- 
and are labels stating that the state has charac¬ 
ter dominated by radial excitation quanta and 
band excitation quanta. Similarly, eigenfunctions of the 
double-well tweezer potential Eq. Q may be written as 
'ipn.p, where p is a vector of parities and n a vector of 
excitation quanta labels. 
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B. Hubbard parameters 

The transition from few- to many-particle physics in 
the presence of a trapping potential is frequently done 
by means of a Hubbard-type model, which projects the 
full many-body model onto a basis of low-energy lattice 
states. Such a projection often removes irrelevant degrees 
of freedom without significantly modifying the physical 
behavior, resulting in models which are easier to analyze. 
The many-body Hamiltonian in second quantization is 


teractions can be well-modeled by zero-range pseudopo¬ 
tentials: an s-wave pseudopotential for bosons 

.Hs-wave (r) = (>') 9rr (11) 

and a p-wave pseudopotential for identical fermions [15] 
Hp-wave (r) = Vd (r) ^ rdrrrr'^ ■ (12) 


^ — j dr‘ip^ (r) 




^(r) 


( 7 ) 


-h 


drdr'ip^ (r) (r') Hi^t (r — rO V' (''O V' (i’ 


where the first line represents single-particle physics in 
the trapping potential V (r) and the second line is two- 
body interactions with interaction Hamiltonian iJjnt (r). 
This Hamiltonian can be expressed in terms of a partic¬ 
ular single-particle basis {ipi, (r)} by expanding the field 
operators in terms of this basis, V'(r) = 
where dj/ is an operator which destroys a particle in 
state ly. A Hubbard model results when this expansion 
is not complete, but runs only over low-energy states in 
the single-particle basis. These basis functions are usu¬ 
ally taken to be localized Wannier functions {wi^ (r)}, 
where i denotes the lattice site (or potential minimum, 
in the case of a double-well optical tweezer) where the 
function is centered, and /i are any other single-particle 
quantum numbers or labels. The precise construction of 
Wannier functions we use will be discussed in Sec. El 
A primary reason to choose the Wannier orbitals as the 
basis is that they decay exponentially in space and the 
couplings can therefore be truncated, for example includ¬ 
ing only nearest-neighbor tunneling in the tight-binding 
approximation. The single-particle contributions to the 
Hubbard model with these approximations can be writ¬ 
ten as 

Hi = Efihifi — + h.c.], (8) 

M (ip) 


Here, Ug is the s-wave scattering length, the p-wave 
scattering volume, and the arrows on the V operators 
denote the direction of operation. Keeping only inter¬ 
actions between particles on the same lattice site, the 
interaction contribution to the Hubbard model can be 
written 








\ M 


U, 




1 + 


Vd^ ,d^, 




_ P-V ' /. 


(13) 


where the s- and p-wave Wannier integrals are 


U, 




= J drw*^,^ (r) w*,, (r) (r) (r) 


(14) 


— dr 


V, (r) (r)) - (r)) (r) 

• (r)) (r) - (r) (r))] . (15) 

Thus, determining the Wannier functions defines an ef¬ 
fective Hubbard model H = Hi + H 2 with Hi and H 2 


given by Eqs. and (13). 


III. NUMERICAL METHODS 


where p. runs over the restricted subset of single-particle 
states, i runs over all lattice sites, (*, j) denotes a sum 
over neighboring lattice sites i and j, dl^ creates a par¬ 
ticle in state p at lattice site i, and 


E, 

J, 


J drw*^{r) 

- J drwp^ (r) 






w, 






(r) 


( 9 ) 

( 10 ) 


with i and j nearest neighbors, are the on-site energies 
and tunneling matrix elements, respectively. 

At ultracold temperatures, only the lowest partial 
waves contribute to scattering. For neutral atoms, in¬ 


In this section we discuss the numerical methods we use 
to obtain the Wannier functions and, from these, the ef¬ 
fective Hubbard parameters for particles trapped in non- 
separable potentials. In particular, we advocate the use 
of discrete variable representations (DVRs) as a flexible, 
simple, and efficient means of solving the single-particle 
problem in anharmonic optical potentials. This choice 
was motivated by the desire to have the very rapid con¬ 
vergence of a spectral method while still maintaining the 
flexibility and simplicity of grid-based methods. Rapid 
convergence is desired both to reduce the computational 
demand of solving fully three-dimensional problems and 
also so that possibly a large number of states can be ac¬ 
curately converged with modest resources. 
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A. Discrete variable representations 

In this section, we briefly review the theory of the two 
DVRs we employ in this work: the sine DVR and the 
Bessel DVR. While the idea of DVRs is quite old, it was 
not until their apparent rediscovery in the 1980s that 
they found broad applicability in chemical and molecu¬ 
lar physics (see e.g. [T^ for a review of this history). In 
spite of their widespread use in chemical physics, DVRs 
have received less attention in the ultracold gases com¬ 
munity (however, see [iTl HB]). Perhaps the greatest 
advantage of DVR methods compared to pure spectral 
methods is their simplicity. As will be shown below, ap¬ 
plying DVR methods requires only the diagonalization 
of a matrix whose elements are all analytically known, 
as contrasted with spectral methods in which the Hamil¬ 
tonian matrix elements consist of integrals which either 
have to be performed analytically for each problem in¬ 
stance or approximated numerically. The advantage of 
DVR methods over other grid-based methods, such as 
finite-order finite differencing schemes, is in efficiency. 
Standard finite-differencing schemes have an error which 
scales as a power of the number of grid points, while DVR 
methods converge exponentially in this parameter. 

We define a DVR as follows m- Given a domain 
I? C K and a Hilbert space H of functions on T>, we 
would like to find a subspace Sk of H with refinement 
parameter K which captures the part of H spanned by 
low-energy eigenstates of our potential. Let Vk be a 
Hermitian, idempotent projection operator into the sub¬ 
space Sk, and pick some set of N grid points {xn}, where 
N = dim5if[3^. We then call the set of Vk and {x„} 
a discrete variable representation if the projected delta 
functions \xn) = Vk^ (x — Xn) are orthogonal. A nor¬ 
malized set of these projected delta functions, which we 
will denote in Dirac notation as |A„) = VK\xn)/VNn, 
with Nn a normalization factor, becomes the set of basis 
functions in which we expand our wave functions of in¬ 
terest. Methods which use an expansion in terms of basis 
functions defined on a grid in real space are also known 
as collocation-spectral methods or pseudospectral meth¬ 
ods [20]. DVRs combine several nice properties of grid- 
based and spectral methods, such as exponential conver¬ 
gence in N of eigenenergies and eigenvectors for problems 
with potentials V (x) G Sk, a diagonal representation for 
the potential energy (A„|V (x) |A„') « V (x„) dn,n', and 
an analytic representation of the kinetic energy operator. 


1. Sine DVR 


1 



x/Ax 

FIG. 3: Sine DVR basis funetions The sine DVR basis 

functions (x|Ao) (red solid) and (x|Ai) (blue dashed) are 
nonzero at their centering grid point and vanish at all other 
grid points. The black dots denote the equally spaced grid 
points. 


where {x\k) = are plane waves with a delta- 

function normalization, and N equally spaced grid points 
Xn = tl'k/K. The projection Vk can also be viewed as 
projection onto the subspace of the Hilbert space spanned 
by non-interacting states with energy less than Ek = 
h?K'^/2M, M being the mass. It can be verified by direct 
substitution that delta functions projected into Sk are 
in fact sine functions 


(a:|A„) 


, sine (tt (x 

TAx 


Xn) /Ax) , 


(17) 


where Ax = Tr/A is the grid spacing, and that these 
functions are nonzero at a single grid point and vanish 
on all others, thus satisfying the orthonormality prop¬ 
erty of a DVR, see Fig. As will be discussed further in 
the section of convergence, the momentum-space cutoff 
K is related to the real-space cutoff Ax, and so the con¬ 
vergence behavior of DVR methods can be interpreted 
either in real or momentum space. 

One notes from Eq. ( [T7| ) that the sine DVR functions 
define a quadrature rule [5^ with abscissas {x„} and 
uniform weights Wn = Ax such that the overlap between 
two such states is exact: 


J dx{An\x) (x|A„/) 





n)) sine (tt (i 


(18) 


^ )) — dn^n' ■ 


The sine DVR |HT| is obtained by letting the domain be 
the real line V = (— 00 , 00 ), choosing the projection op¬ 
erator Vk to project into the subspace of wave functions 
whose bandwidth is limited by a momentum K, 

Vk= f dk\k){k\, (16) 

J-K 


This remarkable property is also the underpinning for the 
Nyquist-Shannon sampling theorem, which states that 
any band-limited function can be completely determined 
from a sequence of equally spaced samples of the func¬ 
tion |22|. Noting that tt/Ax is the spatial Nyquist fre¬ 
quency for a function with bandwidth K, the represen¬ 
tation of a band-limited function in DVR basis functions 
is nothing but Shannon’s interpolation formula. 
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More useful than the fact that we can integrate prod¬ 
ucts of basis functions (and hence, products of any two 
functions in Sk) exactly with the given quadrature is the 
fact that integration by this quadrature produces expo¬ 
nentially accurate results for potentials V (x) which are 
smooth and slowly varying (i.e. well-approximated within 
Sx |40j 1. Here, exponential accuracy refers to conver¬ 
gence with N. The exponential accuracy of the DVR 
quadrature leads to the matrix elements of the potential 
within the DVR basis states, (A„|V (a:) |A„/) = Vnn' ~ 
V (xn) Snn'- As another example of the use of the DVR 
quadrature, the derivative of a band-limited function is 
another band-limited function of the same bandwidth, 
and so the quadrature above provides an exponentially 
accurate representation of differential operators acting on 
functions in Sk via 


(An I 


Qk 

dx^ 


|A„,) 


—-psinc (ttx) 
dx’^ ^ ’ 



(19) 


For example, this gives that the derivative of a function 
expressed in the DVR basis 


^ V>n(a^|A„) 

n 


is given by 


dtp (x) 
dx 


n l^n 


n — I 


tpi{x\An) , 


( 20 ) 


( 21 ) 


and that a representation of the kinetic energy operator 
in the DVR basis is given by 




(An I 




|A„,) 


f ^, n = n' 
2MAx^ I 2 otherwise 

k {n—n'Y ’ 


( 22 ) 


One can show |2T] that these matrix elements also cor¬ 
respond to the limit of an order finite-difference ap¬ 
proximation to the second derivative as A —> oo. This 
provides a useful alternative viewpoint of the sine DVR 
as a limiting case of grid-based finite difference meth¬ 
ods. While in one dimension the resulting Hamiltonian 
matrix is no longer sparse due to the fact that the ki¬ 
netic energy operator is long-ranged in the grid space, 
in higher (tensor product space) dimensions the Hamil¬ 
tonian is sparse in the sense that it is block diagonal 
along each dimension. In addition, while a finite-order 
finite differencing scheme would produce a more sparse 
representation of the kinetic energy than the DVR, such 
methods are only polynomially convergent. Hence, the 
loss in sparsity of the DVR compared to finite-order fi¬ 
nite differencing methods is more than compensated by 
the increase in accuracy. 

A particularly nice corollary of the fact that expec¬ 
tation values of functions well-approximated within Sk 


are exponentially accurate, which seems not to have been 
recognized yet in the literature, is that the DVR repre¬ 
sentation gives an exponentially convergent method for 
the evaluation of integrals of products of eigenfunctions 
and possibly also their derivatives. Such integrals ap¬ 
pear in the evaluation of pseudopotential matrix elements 
Eq. (141-( 15), and their proper evaluation is key to quan¬ 


titative connections between few and many-body physics. 
A demonstration of the convergence of interaction matrix 
elements is given in Sec. IHIA3I 

For potentials with mirror symmetry, V {—x) = V (x), 
we can divide the space of wave functions into those with 
even or odd parity about x = 0. Correspondingly, we can 
use the parity-adapted sine DVR basis sets 


|A+) = 


y/2 (1 -|- Sn,-n) 


(|A„) -|- |A_„)) , n S [0, N] 


(23) 


|A-) = ^ (|A„) - |A_„)) , n e [1, A] , (24) 

which form complete bases for the even and odd functions 
in Sk, respectively. The kinetic matrix elements in these 
bases are 




T+ , = 

2MAx2 I 3 


r ^ _ 2(-i)"+ 

^ {n -D-n ’ 


T~ , = 




2MAx2 1 2(-i) 


{n+n')‘^ 

2(-l)"+"' 
(n-t-nq^ ’ 


(25) 


{n—n'Y 


(n+n') 


n = n 

— , Tl ^ ft 


(26) 


and the potential matrix elements remain unchanged, i.e. 


VI,, = V„, 


V (Xn) S„ 


2. Bessel DVR 


The Bessel DVR is similar to the above sine DVR, but 
uses the free particle wave functions relevant for a radial 
coordinate. Namely, we consider the functions (pm {f) = 
y/rtpm {r) in 2D, where tpm {r) are the solutions of the 
radial Schrodinger equation for angular momentum m € 
Z. The radial functions (pm {f) satisfy the equation 


r d'^(p (r) 
TM dr^ 


rrV — 1/4 


(p{r) 


■VV {r)(p (r) = E(p (r) , 
(27) 


and so behave as near the origin. For free par¬ 

ticles, V (r) = 0, the solutions are Bessel functions: 
(pmk (f) = V^Jm^kr) = {r\km), where k denotes that 
this is the solution with energy lVk^/2m. As with the 
plane wave states, these free radial states obey a delta- 
function normalization. We now obtain a Bessel DVR 
on the radial domain T) = [0, oo) by choosing the projec¬ 
tor VKm to project into the space of free wave functions 
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rjK 


The DVR quadrature associated with the Bessel DVR 
is less useful than that associated with the sine DVR due 
to the fact that the grid points are different for different 
values of the angular momentum m. However, overlaps 
between any two functions which have the same value 
of m still can be evaluated with exponential accuracy 
by this means. The evaluation of interaction matrix el¬ 
ements for radially symmetric potentials is discussed in 
more detail in Appendix [A] 

3. Convergence 


FIG. 4: Bessel DVR basis functions The Bessel DVR basis 
functions (r|A 3 ,i) (red solid) and {r|A 3 , 2 ) (blue dashed) are 
nonzero at their centering grid point and vanish at all other 
grid points. Here, the grid (black points) is not uniformly 
spaced. 


with energy less than Ek = /2M which vanish as 

^m,-i-i /2 origin, 

'PKm = [ dk\km){km\, (28) 

Jo 

and letting the grid points correspond to the zeroes of the 
function (r\Km), the free particle wave function evalu¬ 
ated at the momentum cutoff K. We note that this grid 
is not uniformly spaced, and depends upon the angular 
momentum m. Denoting the n**' zero of Jm (x) as Zmn 
(n = 1,..., oo), the DVR basis functions are 

(r|A„,„) = (-1)” (Kr) , (29) 

and the grid points are ZmniK- Examples of these func¬ 
tions for TO = 3 are shown in Fig. 

As was the case with the sine DVR functions, the 
Bessel DVR functions also define a quadrature which is 
exponentially convergent for functions in Sk- One may 
be concerned about the accuracy of the kinetic energy 
operator evaluated with DVR quadrature, as the cen¬ 
trifugal potential is singular and hence not within Sk- 
However, our DVR basis functions are constructed out 
of the eigenfunctions of the kinetic energy, and so the 
DVR basis also represents the full kinetic energy operator 
with exponential accuracy. Stated differently, the singu¬ 
lar contributions from the derivative operators and the 
centrifugal potential cancel within our basis, and the re¬ 
mainder is well-represented in Sk- The matrix elements 
are given as [23] 

/A I 

( 

_ j 

2M 1 


dr'^ 


-1/4' 




1 11 , 1 . _ 


(- 1 ) 


1 + 

n — n' SZmnZ-rnn' 


n = n 
, n' 


(30) 


In both the sine DVR and Bessel DVR there exist two 
convergence parameters. The first is the number of grid 
points V, and the second is the finite domain size of 
the DVR grid, which can either be taken as a cutoff in 
real space or momentum space. While the domain of the 
DVR grid points is finite, we stress that the DVR basis 
functions themselves still exist on infinite or half-infinite 
domains. The relation between the real space cutoff R 
and the momentum space cutoff K \b R = Ntt/K ior the 
sine DVR and R = ZmN/K for the Bessel DVR. In this 
section, we demonstrate the exponential convergence of 
the sine DVR method with N and R using a ID Gaus¬ 
sian well. The convergence behavior of the Bessel DVR 
method is similar, and we do not demonstrate it here. 

In Fig. [^we demonstrate the exponential convergence 
of the energies for a ID Gaussian potential well 

V(x) = -Ve-'^^^l< , (31) 


where we take V = 96 kHz [41] . Wo = 707 nm, and 
the mass of ®^Rb. This convergence is monitored by in¬ 
creasing a convergence parameter r\ in discrete steps d, 
and then observing an exponential decrease in AE(r]) = 
E(r] + 6) — E(r]). Panel (a) demonstrates exponential 
convergence with the number of DVR grid points N at 
fixed domain size R, which can also be stated as conver¬ 
gence in the grid spacing Ax. The different curves denote 
different eigenstates, with lower curves corresponding to 
lower energy eigenstates. In panel (b) we show exponen¬ 
tial convergence in the domain size R at fixed Ax. In this 
case, convergence of the hrst 10 eigenenergies to machine 
precision can be achieved with only 60 grid points (in the 
parity-adapted DVR basis set). 

In Fig. we demonstrate the convergence of the one¬ 
dimensional s- and p-wave interaction integrals, Eq. (14) 
and Eq. ( jl^ , for the same parameters as Fig.[^ Namely, 
we plot the differences in the dimensionless parameters 
Uo,o,o,oo-ho, C7io,io,io,ioaho, Goqo.io.oOho, and Vbqouo.oOho 


with neighboring convergence parameters, where Oho = 
^/h/rnuj is the harmonic oscillator length corresponding 
to the trap curvature. As with the energies, we use differ¬ 
ences between neighboring convergence parameters, e.g. 

= (Uofififi(N + SN) — Uofififi(N))aYio to 
gauge convergence. We have also checked the conver¬ 
gence of the interaction matrix elements for the harmonic 
oscillator, where analytic results are available, and found 















Aa;(wo) 



N 




FIG. 6: Convergence of s- and p-wave interaction matrix 
elements in the DVR basis size. The convergence of the one¬ 
dimensional analogs of Eq. (141 and Eq. (151 in oscillator 
units is given by considering the differences with neighboring 
convergence parameters. These integrals show a similar rate 
of convergence to the energies in Fig. [^a), though the expo¬ 
nential convergence sets in at a larger value of N. The trap 
parameters are the same as those used in Fig. 


B. Optimizations: Using Sparsity and Variational 
and Quasi-adiabatic methods 


FIG. 5: Convergence of DVR energies with number of DVR 
grid points and domain size, (a) Gonvergence in the num¬ 
ber of DVR grid points (equivalently, the grid spacing A®) is 
demonstrated for the first 10 eigenstates of a Gaussian well 
at fixed R = 3wq. The points are differences in energies for 
neighboring N, as described in the main text. Curves from 
bottom to top are increasing in energy, (b) The analogous 
plot to (a) at fixed Aa; = 0.05 and varying R, the DVR grid 
domain cutoff. The parameters used are: depth V = 96 kHz, 
waist Wo = 707 nm, and the mass of ®^Rb. 


similar convergence. The rate of convergence of the ma¬ 
trix elements is akin to that of the energies. However, 
exponential convergence sets in at a larger value of N for 
interactions compared to the convergence of the energy. 
Convergence in R has a similar qualitative behavior. Due 
to the non-uniform grid of the Bessel DVR, computation 
of the interaction parameters for radial functions is more 
involved, and is discussed in Appendix [X} 


In order to converge results, we find it is useful to first 
converge in Aa; and then increase R until convergence. 
Convergence with Ax can be ascertained by convergence 
of the energy at fixed R, and then convergence with R can 
be judged by requiring that the weight on the outermost 
DVR basis functions becomes on the order of machine 
precision. In higher-dimensional scenarios, visualization 
of the convergence with R can be greatly assisted by plot¬ 
ting the nearest separable state. A precise definition of 
the nearest separable state and a procedure for calculat¬ 


ing it are discussed in Sec. IV 


A nice feature of DVRs in dimensions higher than one 
are that the resulting multi-dimensional Hamiltonian de¬ 
scriptions are sparse. In particular, for a 3D problem in 
Cartesian coordinates, the Hamiltonian may be written 
in a product sine DVR basis as 

(a,,Aj,a^|h|a;a;a;) = 

~V ^x.x'^y^y’^zz' ~V ^x,x'^y ,y'^ z.z'^ 

Hence, the Hamiltonian can be applied to a vec¬ 
tor describing the state in the DVR basis us¬ 
ing O {^N'^NyNz + NxNyNz + NxNyNf) operations, far 
fewer than the O [N^NyN^) operations required for a 
dense matrix description. This enables the use of sparse 
matrix methods requiring only matrix-vector products, 
such as the Lanezos algorithm [21], to find extremal 
eigenstates. Utilizing the sparsity of the representation 
is also key for efficient simulation of the dynamics of par¬ 
ticles in time-varying potentials, as it enables the use of 
Krylov subspace approximations to the matrix exponen¬ 
tial (23 US]. Such dynamical simulations are useful, for 
example, for determining the degree of adiabaticity when 
optical tweezer wells are dynamically repositioned. 

For lattice problems, using a set of basis functions 
which has the same translational symmetry as the lat¬ 
tice often leads to significant computational gains. The 
DVR basis functions given above are defined on infinite or 
semi-infinite spaces, and so are inappropriate for expand¬ 
ing a function defined on a finite, periodic space. For ID 
periodic potentials, expansion in terms of a plane wave 
basis is very efficient and accurate. Hence, it is natural 
to combine DVR methods for the transverse potential 
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with plane-wave methods along the periodic direction. In 
what follows, we combine the two by expanding the full 
Hamiltonian in a basis consisting of products of trans¬ 
verse (r) and lattice (z) degrees of freedom. The number 
of functions we use for the expansion along, e.g. , direc¬ 
tion r will be denoted a^, and is called the variational 
dimension. 

To see how the combination of DVR methods with 
other basis sets is facilitated in this situation, consider 
the lattice potential Eq. ([^. We can write this poten¬ 
tial as EiattVg (r) 14 (z), where Vg (r) = -exp(-^) is 

the Gaussian radial potential and 14 {z) = cos^ (kz) is 
the lattice corrugation. This multiplicative separability 
of the potential provides a natural variational basis with 
which to expand the full coupled 3D problem, namely 
products of eigenfunctions of ViattVg (^) with eigenfunc¬ 
tions of VuttVz (z). Let us denote the eigenfunctions 
of Vg{r) with angular momentum m as 'R-m,nr. (’") and 
the eigenfunctions of Vz{z) with quasimomentum q as 
(z), where and Uz are mode and band indices. 
Then, the matrix elements of the full Hamiltonian in the 
product basis (rln^n^) = (r) {z) 

are 


+ (33) 


1(2) 


where is the radial and azimuthal 

kinetic energy operator, and is the z kinetic energy 
operator. If we use the Oj. ^ N.^. lowest-energy states 
\nr) to expand the full coupled problem, where we will 
call Or the variational dimension and Nr is the DVR ba¬ 
sis size, can be efficiently obtained with O (a^V^) 

operations within a DVR calculation, as the matrix V is 
diagonal in the DVR representation. The kinetic energy 
matrix elements can then be obtained with O oper¬ 
ations as Tnn' = dnn'En-Vnn', being the eigenenergy 
of state n. Restricting to a variational basis results in an 
eigenvalue problem of linear dimension Orttz whose low¬ 
est energy solutions converge to the true solutions with 
Oryttz —> 00 . The advantage of this procedure is that 
the variational dimensions Or and can be significantly 
smaller than the “bare” basis sizes Nr and Nz required 
to converge the variational basis functions TZ and ip. This 
procedure can also be applied to the potential Eq. §. 
Here, we take the basis states TZ to be the eigenvectors of 
the Hamiltonian with potential {I4onst + V) V® (r) and 
the states (p to be the eigenvectors of the Hamiltonian 
with potential VV^ (z). The variational Hamiltonian has 
matrix elements 


{nrnJ\H\nrnz) — E 

E {Vconstdn^n'^ + ■ ( 34 ) 


In principle, one could optimize the depths of the po¬ 
tentials used to compute the functions TZ and p such 
that the variational dimensions and are minimal, 



FIG. 7: Convergence of the m = 0 energies and tunneling 
amplitudes of the potential Eq. § with the variational basis 
size Oz. The difference in energies adding another lattice band 
AE = E{az + 1) — E{az) as a function of eigenstate energy are 
shown for numbers of bands Oz = 2,..., 7, with lower curves 
corresponding to larger numbers of bands. The bottom curve 
is AJ for Oz = 7, showing that the tunneling amplitudes are 
display similar convergence behavior to the energy. The trap 
parameters are wo = 30 pm, Konst = 88 Er, and V = 12 Er. 


but in practice using the depths given above works well 
and we have not explored this possibility. We demon¬ 
strate the convergence of this variational procedure in 
the number of bands kept with the potential Eq. ([^ 
with wq = 30 pm, Konst = 88 Eh, V = 12 Ep, where 
Er = /2ma'^ is the recoil energy of ®^Sr in a magic 

lattice of spacing a = 406.72 nm. Fig. displays the 
results, showing rapid convergence in both the energies 
and tunneling amplitudes. The transverse variational di¬ 
mension is oo 500. For comparison, the DVR basis size is 
~ 2000, and the number of plane waves used to expand 
the lattice potential is ~ 200. 

For optical tweezers, in which the Rayleigh range can¬ 
not be neglected, the potential no longer becomes a a 
product of functions in the transverse and z degrees of 
freedom. In this case, one can employ a quasi-adiabatic 
variational method, in which a set of basis functions are 
constructed for each axial (z) DVR basis state |Az^) us¬ 
ing only the transverse kinetic energy, and then the full 
problem is diagonalized in this basis. As an example, 
consider the single-well optical tweezer potential Eq. (HI- 
We construct the basis functions \nr, Zn) as the Or lowest 
energy eigenfunctions of T^^'> + K {r,Zn) for each DVR 
grid point z„, and denote the eigenenergies as En.,.^z„- 
These energies then form an effective potential for the z 
degrees of freedom, which obey the Schrodinger equation 

(n,., Zyi |n^, Zy,) - E En.,.^Zri^Ur,n'.^^Zn,z!^ ■ 

(35) 

The time required for this complete procedure scales as 
O {NzD{Nr) -b D {urNz)), where D{Ng) is the time re¬ 
quired to find the lowest eigenvectors of a (possi¬ 
bly sparse) matrix of linear dimension N^. As with the 
above, reductions in the computational time required for 






10 


a given accuracy can be achieved when ar ^ Nr. 

A similar procedure can also be applied to the double¬ 
well optical tweezer potential Eq. Q, using the multi¬ 
plicative separability of the x and y potentials. Here, 
we first diagonalize the x kinetic energy together with 
the potential Vd (xjO, z„) at each axial DVR grid point 
Zn, keeping the lowest quasi-adiabatic levels. Next, 
we diagonalize the x and y kinetic energies together 
with the full potential constructed from the product 
of the matrix elements of Vd {x, 0, Zn) in the x quasi- 
adiabatic basis and a Gaussian along y, extracting the 
lowest ttxy quasi-adiabatic energies. These energies then 
form an effective potential which is diagonalized with 
the 0 kinetic energy to obtain the full solutions in 
O {N^D{Nx) + N^D{axNy) + D{axyN^)) time. 


IV. MATRIX PRODUCT STATE 
REPRESENTATIONS OF NON-SEPARABLE 
STATES 


The eigenfunctions of a non-separable potential are 
themselves not separable, meaning that ^ (r) ^ 

tpx {x) ipy {y) tpz {z) or V' (r) ^ {(t>) A (r) tpz (z) for po¬ 

tentials of cylindrical symmetry. However, in many cases 
one would expect that the lowest-energy states of sufh- 
ciently deep potentials are “nearly separable,” as the po¬ 
tential becomes nearly harmonic over the extent of the 
wavefunction. In this section, we provide several quanti¬ 
tative measures of non-separability for the eigenstates of 
non-separable potentials, and investigate near-separable 
approximations of non-separable states using tools bor¬ 
rowed from quantum information theory. 


A. Schmidt form for non-separable states with 
cylindrical symmetry 

The simplest case to study non-separability in a 3D 
potential is to consider a system with cylindrical symme¬ 
try. Here, the azimuthal degrees of freedom are separa¬ 
ble, leaving only the r and z degrees of freedom coupled 
so that the wave function may be written 

V'm (r) = -^Cm {r, z) . (36) 

We can find the state nearest to the true state in the 2- 
norm with a restricted amount of separability (in a sense 
to be made precise below) by using the Schmidt decom¬ 
position m- The Schmidt decomposition states that we 
can write any state on a product Hilbert space T-Li® ^2 
as 

W = (37) 

fi=i 

where the sets {\(j)fi)} and {|Cm)} consist of orthonormal 
states in "Hi and 7^2, respectively, and the Schmidt co¬ 
efficients {A^} satisfy > 0, = 1- Schmidt 


decomposition is unique, up to unitary rotations in sub¬ 
spaces with degenerate Schmidt coefhcients. The dimen¬ 
sion X is called the Schmidt rank, and % = 1 if and only 
if the state is separable. Hence, x quantifies the degree of 
non-separability. A continuous measure of separability is 
given by the von Neumann entropy S = — A^ log A^. 
We note that the Schmidt decomposition can be obtained 
efficiently numerically using the singular value decompo¬ 
sition (SVD) of the coefficient matrix of the wavefunction 
in a given basis, e.g., a DVR basis. Further, the Schmidt 
decomposition is invariant under local unitary transfor¬ 
mations, which is to say that Eq. (37) is unchanged if we 


change bases in or 772- This implies that if we have 
our wavefunction expressed in terms of a variational basis 
formed of products of r and z functions, as described in 
Sec. |HIB[ then we can find the Schmidt decomposition 
by performing the SVD on the matrix of coefficients in 
the variational basis. As the variational basis sizes a^, 
are often significantly smaller than the grid sizes Nr, N^, 
the resulting SVD computation is more efficient. 

The state nearest to \^p) with fixed non-separability, 
i.e. fixed Schmidt rank x ^ X: is obtained by truncating 
the expansion Eq. (37) at /i = x- The difference between 


this truncated state and the true state is given by the 
discarded Schmidt coefficients 


= E 

Q:=x+1 


•X • 


(38) 


In practice, the state |'0) needs to be renormalized by 


setting the Schmidt coefficients A^ = \jY^^=i^'on 

H = 1,..., X- The resulting correction of (1 — e,^) 


-1 


IS 


inconsequential to the error bound of Eq. (38) as ev —)■ 0. 


For states with a low degree of non-separability such that 
the Schmidt rank required to reproduce the state with er¬ 
ror e in the sense of Eq. (38), x, is much less than the 


maximum allowed Xj using the Schmidt form drastically 
reduces the memory required for storing quantum states 
from NrNz, to x {Nr + A^z), and also reduces operations 
acting only on one degree of freedom, say the z degree 
of freedom, from NrN^ to X-A^z- The savings are es¬ 
pecially beneficial in cases where thermal averages over 
a large number of states are to be performed. Finally, 
we note that the nearest separable state is obtained by 
setting X = 1- In addition to being useful for visualiza¬ 
tion and developing intuition, the nearest separable state 
IV^sep) dehnes an additional measure of non-separability 
£ = — log(|(V'sep|V')P) which we will call the geometric 
non-separability, in analogy with the geometric entangle¬ 
ment [28] motivating its definition. As opposed to the von 
Neumann entropy, which is a measure of non-separability 
for a specific bipartition of the degrees of freedom, the ge¬ 
ometric non-separability is a global measure of the non¬ 
separability of the full state. This distinction is most 
important in the multi-partite case discussed in the next 
section, where multiple possible bipartitions exist. 
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B. Matrix product state form for non-separable 
states in Cartesian coordinates 


The above analysis is directly applicable in the case of 
a single-well optical tweezer or a ID optical lattice, but 
not for a system without cylindrical symmetry, such as 
the double-well optical tweezer Eq. Q . Here, the state is 
most directly represented in Cartesian space, and the a:, 
y, and 2 degrees of freedom are all coupled. The natural 
generalization of the Schmidt decomposition to multipar¬ 
tite systems is a matrix product state (MPS) decomposi¬ 
tion [T^ (also called a tensor train decomposition in the 
mathematics literature m), which in the present case 
reads 

Xxy Xyz 

V' (r) = XI X (y) (^) ' (39) 

Pxy —1 


Here, Xxy and Xyz are the Schmidt ranks corresponding 
to bipartitions of the state into x and y^z degrees of free¬ 
dom and x^y and z degrees of freedom, respectively. The 
savings in using the MPS format for a non-separable state 
is even more striking than in the two-partite case: storage 
is reduced from Ny;NyN^ to iXxyNx+XxyXyzNy+XyzN^) 
and operations on, e.g., the x degree of freedom can be 
applied in N^Xxy operations rather than N^NyNz- We 
note that the representation Eq. (39) is not the only MPS 
topology which can represent a non-separable state in 3D. 
Por example, we could have chosen the partition x — y — z 
or z — x — y. These other permutations will generally have 
different Schmidt ranks y, which implies a difference in 
accuracy for describing the state for a given number of 
parameters. Here, we do not make any claims about the 
optimality of the x — y — z partition of degrees of freedom 
in Eq. (39); the optimal partition of degrees of freedom 
must be determined on a case-by-case basis. 

The multi-partite representation Eq. (391 bears many 
similarities with the two-partite Schmidt decomposition 
Eq. ( |37[ ). Eor example, both approximations are con¬ 
trolled by the Schmidt ranks of bipartitions, and both 
representations can be obtained via the SVD (recur¬ 
sively applied at each bipartition, in the case of the MPS 
form), as is shown explicitly in Appendix [b| However, 
there are key differences between the MPS form and 
the Schmidt form. An important example is in find¬ 
ing the nearest state with a restricted amount of non¬ 
separability. Compression of an MPS to the optimal MPS 
with restricted non-separability y at a given bipartition 
is achieved by trun cating the range of y,yz in Eq. (B2) 
or pxy in Eq. (B4). The resulting error in the 2-norm 
is given by the sum of squares of the discarded singular 
values, as in Eq. (38). However, finding the nearest MPS 
in a global sense is not as simple as for the two-partite 
case [soiin]. Here, variational algorithms which mini¬ 
mize the distance between the true state and one in a 
variational manifold with fixed non-separability in an it¬ 
erative fashion perform well. In Appendixwe present 
a such a variational algorithm for hnding the product 


state X{x)Y(y)Z{z) nearest a given state expressed in 
the MPS format. 


While the above decompositions are useful for visual¬ 
ization, quantification of non-separability, and storage of 
states, obtaining such a decomposition via the SVD is 
as computationally demanding as finding the eigenstates 
themselves. However, the fact that the states explored in 
this work are only weakly non-separable in spite of the 
fact that they are significantly anharmonic hints that the 
MPS form given by Eq. (39) could be useful as an ansatz 
which is variationally optimized, similar to the way MPSs 
are used in the density matrix renormalization group 
(DMRG) method of condensed matter physics [32] . Such 
an algorithm would find states with a restricted degree 
of non-separabaility using significantly fewer resources 
than direct diagonalization. This prospect is especially 
promising for multi-particle systems, where failure of the 
center of mass and relative coordinates to separate in an 
anharmonic potential makes analysis significantly more 
difficult |33J [31] . We leave a detailed description of such 
an algorithm for future work. 


V. APPLICATION TO OPTICAL TWEEZER 
ARRAYS AND LATTICES 


In this section, we apply the above algorithms and 
measures to two cases of interest. The first is a double¬ 
well optical tweezer of the form Eq. (|^ , as has been con¬ 
sidered in recent experiments at JILA [7]. Here, we focus 
on the properties of the lowest few states. The second 
case we study is that of a ID optical lattice of the form 
Eq. Here, we are interested in the spectrum of states 
spanning a large fraction of the lattice depth, as is of¬ 
ten the case when non-degenerate gases are trapped in 
such lattices. The reason we study the generalization 
Eq. § over Eq. ([^ , with the latter corresponding to the 
potential used in neutral atom optical clocks [11] [14] , is 
that the potential Eq. (§ enables us to study tunneling 
in a non-separable lattice with a tunable degree of both 
non-separability and transverse confinement. 


A. Optical tweezers 


As an example of our methodologies applied to ar¬ 
rays of optical tweezers, we consider a double-well optical 
tweezer of the form Eq. (ffl, and take parameters similar 
to the JILA experiment |7] , with waist wq = 707 nm, 
Rayleigh range zr = 2.17 /xm, and mass and s-wave 
scattering length for ®’^Rb in the \F = 2,mp = 2) 
state. Before we give results for Hubbard parameters 
and quantitative measures of non-separability, we first 
detail a procedure to obtain localized Wannier-like func¬ 
tions from the numerical eigenstates of the double-well 
optical tweezer. 
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1. Construction of Wannier functions 

For inversion-symmetric lattice potentials in ID, the 
localization properties of Wannier functions can be ana¬ 
lyzed in detail. In particular, a seminal work of Kohn m 
showed that requiring the phases on the Bloch functions 
ipnq (x) to be such that ipnq (x) is a smooth function of 
the quasimomentum q leads to Wannier functions which 
are real and have an asymptotically exponential decay 
away from their centering position. Apart from this situ¬ 
ation, a procedure to obtain so-called maximally localized 
Wannier functions has been extensively developed which 
minimizes the second moment of Wannier function posi¬ 
tion away from its center |dflj . Here, we put forth a sim¬ 
ilar criterion for the generation of localized orbitals from 
the eigenstates of the non-separable double-well optical 
tweezer potential Eq. Q. With a slight abuse of termi¬ 
nology, we will also refer to these localized orbitals as 
Wannier functions. 

In particular, we take as our localization functional the 
probability to measure a particle in the left well. Math¬ 
ematically, we use the functional 

/ P^cut 

dydz / dx (r) (r) , (40) 

J —OO 

where Xcut is some parameterization of what defines the 
left well, e.g. the center between the minima, the local 
maximum, or a; = 0. Now, given some subset of spatially 
overlapping eigenstates {ipu^ (r), V'fea (j"): ■ • ■ : V’feiv (i")} 
we seek the normalized linear combination w (a, r) = 

Y.i=i “iV'/ci (r), a • a = 1 satisfying 

max£ [w (a, r) , u; (a, r)] . (41) 

a 



FIG. 8: Single-particle Hubbard parameters in an asymmet¬ 
ric double-well optical tweezer. The maximal localization pro¬ 
cedure described in the text converts the tunneling doublet 
whose energy splitting (e 2 —£i) shown in magenta into a single¬ 
band Hubbard model with a nearly-constant tunneling J in 
red and an effective bias {E 2 — Ei) which scales linearly with 
the applied potential bias. 

for greater efficiency. In the case of a symmetric double 
well, A = 0, the closely spaced doublets with even and 
odd parity along the tunneling direction form a 2 x 2 lo¬ 
calization matrix with diagonal elements 0.5 and small 
but non-zero off-diagonal elements. Hence, the Wannier 
functions in this case are even and odd superpositions of 
the even and odd parity functions, which are related to 
one another by the parity transformation. 

2. Tunneling and effective bias in an asymmetric 
double-well optical tweezer 


Alternatively, for functions maximally localized on the 
right, we take the minimum. Eq. (411 can be written as 


max a • L • a , 

a 


(42) 


where the localization matrix has matrix elements 


Ly = C [ipk, (r), tpkj (r)] . (43) 

L defines a symmetric positive semidefinite quadratic 
form, and so the normalized vector maximizing the lo¬ 
calization criterion is the eigenvector corresponding to 
the largest eigenvalue. Moreover, the complete set of 
orthonormal eigenvectors of the localization matrix pro¬ 
vides a new basis which is optimal according to our lo¬ 
calization criterion. Strictly speaking, in order to find 
the optimally localized basis according to our criterion, 
all states at experimentally relevant energies should be 
included. However, for states which are well-separated 
in energy, the mixing between them incurred in the lo¬ 
calization transformation is small and has little effect on 
the Hubbard parameters. Hence, one can localize subsets 
of energetically separated states (with energy differences 
large compared to the tunneling splitting) individually 


Using the localization prescription given in the last sec¬ 
tion, we find the tunneling amplitudes and on-site en¬ 
ergies Efj^, cf. Eqs. ([^-(10 1 in the double-well in the basis 
of maximally localized Wannier functions. In particular, 
we have 


E^ = Y,{ai^'>)h., (44) 

Jn = , (45) 


where is the eigenvector of the localization ma¬ 
trix, Cl, is the energy of state and /2 denotes the 

index of the state connected to /x by tunneling. 

As an example. Fig. shows the results for an a = 
853 nm spacing and U = 96 kHz depth double-well 
tweezer as a function of the applied bias A, focusing 
on the lowest two states. These two states are a tun¬ 
neling doublet when A = 0, and asymptotically become 
two localized states with negligible spatial overlap as A 
becomes large. The tunneling J(A) in the optimized A- 
dependent Wannier basis changes by less than a percent 
over the range of biases shown here. The effective bias 
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FIG. 9: Tunneling and interaction parameters in a sym¬ 

metric double-well optical tweezer. The tunneling and s-wave 
interaction energies for the ground (g) and first 2 : excited (e) 
states in a double-well optical tweezer as a function of depth. 
While the reduction in Use compared to Ugg is an effect of the 
different spatial extent of the Wannier orbitals in the g and e 
states that occurs even for separable potentials, the increase 
of Je vs. Jg is solely a manifestation of non-separability. 


Aeff, given by the difference in energy (E 2 — Ei) between 
the Wannier functions localized in the two wells, is well- 
represented by a linear function of the potential bias. 
However, the slope of the linear dependence is generally 
less than one, Aes = mA, m < 1, due to curvature 
effects as the potential depth is changed. In the given 
example, the best fit gives m = 0.79. Also displayed 
is the difference in energy between the two eigenstates 
used to construct the Wannier functions, (e 2 — ei). This 
energy difference fits well to (2J)^ -I- A^^j, validating 
the single-band Hubbard model description, even when 

A > J. 


3. Tunneling and interaction parameters in a symmetric 
double-well optical tweezer 

We now turn to the case of of a symmetric double well 
potential, using as an example the spacing a = 820 nm. 
In Fig. we show the tunneling and interaction ener¬ 
gies for the ground {g) state as well as the first axi¬ 
ally excited state along 2 (e). We note that the exci¬ 
tation in the state e is perpendicular to the tunneling 
direction, and so would not affect tunneling in a sepa¬ 
rable potential. The excited state tunneling is roughly 
20 % larger than the ground state tunneling as a con¬ 
sequence of non-separability. To make a quantitative 
comparison with tunneling in optical lattices, we can ap¬ 
proximate our double-well tweezer as an optical lattice 
with lattice constant a twice the distance from the origin 
to mina;>o V (x, 0,0) and depth given by the local maxi¬ 
mum V = V (0,0,0) — V (a/2,0,0). For the given spac¬ 
ing we consider, a = 642.2 nm, the associated energy 
Ea = h^TT^/2ma^ is 1.39 kHz, and the effective lattice 
depth is F = 0.0473F. We fit our data in the range 


FIG. 10: Non-separability of states in a symmetric double¬ 
well optical tweezer. The von Neumann entropy S of xy and 
yz bipartitions and the geometric non-separability £ for the 
ground state (g) and the first 2 excited state (e). While both 
states are near-separable, the excited state is less separable 
than the ground state. 


V G [50 kHz, 200 kHz] to the standard formula for tun¬ 
neling in an optical lattice 




exp (^-C^V/En) , (46) 


using Ea as the recoil energy and V as the lattice depth. 
In contrast to the optical lattice, in which the optimal 
parameters are A = 1.363, B = 1.057, and C = 2.117, 
we find A = 2.563, B = 1.217, and C = 2.281. It should 
be noted that these values also depend on the waist and 


Gaussian spot spacing. While the function Eq. (46) pro¬ 


vides an excellent fit, the optical lattice analogy itself 
amounts to a 60% — 70% discrepancy over the range of 
parameters considered, with the tweezer having a larger 
tunneling for a given effective lattice depth. 

Fig. [9] also shows the behavior of the s-wave inter¬ 
action energies 11^^' = UaT-Ta as a function of 

lattice depth. For an intuitive understanding of the 
behavior and to quantify the degree of anharmonicity, 
we will compare with the predictions obtained by us¬ 
ing a harmonic approximation for the potential minima. 
For harmonic wells, U ~ (axOi/Oz)”^, with Ojy the har¬ 
monic oscillator length along Cartesian direction v, and 
Qi, ^ leading to U ^ The best fit for 

Ugg predicts Ugg ~ yo.85^ showing a faster rise of in¬ 
teractions with lattice depth than predicted for a har¬ 
monic well. Further, in a harmonic potential, the inter¬ 
actions are related as Uee/Ugg = 3/4 and Ueg/Ugg = 1/2, 
but due to anharmonicity Uee/Ugg « 0.714 — 0.728 and 
Ueg/Ugg « 0.484 — 0.49 for the range of depths in Fig. 


4 . Separability characteristics of double-well tweezer states 

Fig. shows the von Neumann entropy of non¬ 
separability for bipartitions into x and y ® z degrees of 
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FIG. 11: Nearest separable states of a double-well opti¬ 

cal tweezer and their harmonic approximations. The compo¬ 
nents of the nearest separable wavefunction X{x/wo) (red), 
Y{y/wo) (green), and Ziz/zn) (magenta) for a symmetric 
double-well optical tweezer with depth = 52 kHz are 
shown as functions of their respective dimensionless coordi¬ 
nates. The harmonic approximations are also shown for y 
(dark blue) and 2 (light blue), demonstrating significant an- 
harmonicity in the 2 degree of freedom. 


freedom (xy) and x®y and 2 degrees of freedom (yz), as 
well as the geometric non-separability, for the states given 
in Fig.|^ In spite of the fact that properties, e.g. the tun¬ 
neling in the excited state, show significant characteris¬ 
tics of the non-separability of the potential, each eigen¬ 
state itself is very nearly separable, requiring a Schmidt 
rank of y = 5 across any bipartition to capture the state 
in the 2-norm to an accuracy of 10“^^ (see Eq. (381). This 
near-separable character is also borne out in the separa¬ 
bility measures of Fig. (§. For example, the geometric 
non-separability £ demonstrates that the ground and ex¬ 
cited states can be represented by their respective nearest 
separable approximations with fidelities of « 0.9995 and 
« 0.9987. As a general rule, bound states of higher en¬ 
ergy are less separable than bound states of lower energy. 
Also, it is interesting to note that the degree of separa¬ 
bility of the states of the double-well optical tweezer are 
not monotonic functions of the tweezer depth. 


The components of the nearest separable state if (r) = 
X{x)Y{y)Z{z), obtained using the method in Ap¬ 
pendix are shown in Fig. El In addition, the har¬ 
monic oscillator approximations for the y and z degrees of 
freedom, e.g., Fho(y) = exp{-y^/{2alJ )/with 
Oho = wo{2V/Ei^g)~^/'^, are shown for comparison. The 
nearest separable state along the y direction is accurately 
represented by its harmonic approximation, while the 
2 state is considerably wider than its harmonic coun¬ 
terpart and has a different shape, showing strong an- 
harmonicity. Finally, we stress that even though two 
states ip{r) and may both be nearly separable, 

'i/’(r) « X.^{x)Y^(y)Z^{z) and ^(r) « X^{x)Y^{y)Z^{z), 
this does not imply that Ay,(a:) « X^{x) etc.. 


B. Optical lattice 

In this section, we turn our attention to a non- 
separable lattice of the form Eq. © , where we will fix 
V = \2Eji, with En = /2ma? the recoil energy 

of ®’^Sr in a lattice of spacing a = 406.72 nm, and 
Wq = 30 /im. The “running-wave” component of the 
lattice Fconst will be left as a variable to show the effect 
of increasing transverse confinement on the axial motion. 
The transverse confinement can be characterized by the 
effective oscillator frequency fux = (Fconst + F) 

resulting from a harmonic expansion near the Gaussian 
potential minimum. 


1. Construction of Wannier functions 


As discussed in Sec. IIIB the states of the non- 


separable ID optical lattice may be variationally ob¬ 
tained in the form 


(l*) — 






/ , Cn-r 


nn 


zE'm.n (^) ^quz ('^) : ( 47 ) 


where 77.(r) are a set of functions obtained from a DVR 
calculation using the transverse Gaussian potential and 
Lpqriz (-z) are Bloch functions with quasimomentum q and 
band index obtained from a plane-wave calculation. 
Here, we briefly describe how we construct localized 
Wannier functions in this non-separable case by anal¬ 
ogy with the separable case worked out by Kohn [TU] . In 
Kohn’s original scenario, which is an inversion-symmetric 
ID lattice, real Wannier functions of maximal local¬ 
ization are obtained by using the transformation prop¬ 
erty of Bloch functions under inversion 2 ) = 

(—( 2 ) (we index starting from 1) and the 
requirement that the Bloch functions are smooth func¬ 
tions of q to fix the gauge, i.e. phase, ambiguity in the 
Bloch functions. Using this choice of phases, maximally 
localized Wannier functions follow from 


Wnzi{z) 


1 

71 


(jeBZ 


(48) 


where L is the number of lattice sites and Zi the centering 
site location. The resulting Wannier functions transform 
under inversion as Wnz,i{—z) = the 

Wannier function center is inverted and a phase may be 
acquired. 

In the non-separable case, we can generalize the 
inversion symmetry transformation to ipqi, (r, — 2 ) = 

P^ip-qv (r, 2 ), where = 1 if the dominant weight of 
the state Tp lies in bands = 1, 3, 5,... and —1 other¬ 
wise. This phase can also be set unambiguously at the 
BZ center q — 0, where translations and inversions com¬ 
mute. The remaining phase ambiguities are the Bloch 
function phases in q € (0,7r/a), which become ±1 un¬ 
der the requirement of real Wannier functions. We fix 
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FIG. 12: Tunneling vs. Energy. The tunneling, in units of the 
separable lattice tunneling, as a function of Wannier state en¬ 
ergy for V = 12 En and uiq = 30 /rm. The upper pair of curves 
corresponds to Vconst = 88 Er , giving a transverse frequency 
huj ~ 400 Hz, with the red (green) curves corresponding to 
azimuthal quantum number m — 0 (m — 200). The lower 
pair of curves has V);onst = 388 Er {huu ~ 850 Hz), with blue 
(magenta) corresponding to m = 0 (m = 200). J depends 
dominantly on the transverse mode energy, as intuitively un¬ 
derstood through the lowest-order harmonic expansion of the 
transverse potential. 


the gauge here by requiring that the Bloch functions be 
smooth functions of q, by analogy with Kohn’s original 
work. With these phase conventions, Wannier functions 


are constructed from Eq. (471 as 


^^ (jGBZ 




mqu 


(r) . 


(49) 


Our construction produces Wannier functions which 
transform into each other up to phases with the symme¬ 
tries of the lattice and reduce to the maximally localized 
Wannier functions of Kohn when the transverse motion 
separates from the axial motion. We note that the partic¬ 
ular choice of Wannier functions affects only interaction 
matrix elements and not tunneling matrix elements, as 
the latter are set by the energy as a function of q alone, 
see Eq. (pi. 


2. Tunneling and interaction parameters 

Using the energy dispersion Em,q,v of state Z’mqv (r), 
the tunneling of Wannier state Wmvi (r) between neigh¬ 
boring lattice sites is given by 

J™,. = -;^ (SO) 

(jGBZ 

i.e. the Eourier transform of the band structure in the rel¬ 
ative lattice coordinate. Similarly, the Wannier function 
energy, Eq. (§, is given as the average of the dispersion 


in the BZ, 


Em.if — ^ ( EjYi q i, . 


(51) 


ijeBZ 


In Eig.[^we show the behavior of the lowest-band tun¬ 
neling as a function of the Wannier state energy, with 
the former measured in units of the separable (ID) lat¬ 
tice lowest-band tunneling Jgep and the latter measured 
in units of the recoil energy, with the zero of energy being 
the ground state energy. We define the lowest band as 
being the set of states such that Eq. (47) 

is maximal for = 1, where riz labels the “bare” bands 


used in the variational expansion of Eq. (47). For the 


parameters we consider, the mixing between the bare 
bands is slight (> 85% of the population is in = 1 
for the energy range we consider) and there is no ambi¬ 
guity in this definition. Hence, the energy on the cc-axis 
of Fig. correlates dominantly to transverse mode en¬ 
ergy. The tunneling generally increases with increasing 
energy, which can be understood by expanding Eq. (§ 
to lowest order in r 


%latt — Ujonst ( — 1 + 2 


Wr 


V \-l 


wt 


cos^ {kz) . (52) 


Hence, to lowest order, the lattice potential depth is low¬ 
ered by an amount proportional to (r^), which is itself 
proportional the transverse mode energy in the harmonic 
oscillator approximation. Semiclassically, a particle in a 
higher transverse mode spends more time near the classi¬ 
cal turning points, and here the potential depth is smaller 
due to the non-separability. Because of the non-linear re¬ 
lationship between lattice depth and tunneling, the rela¬ 
tionship between transverse mode energy and tunneling 
is also non-linear. The semiclassical reasoning for the 
dependence of tunneling on energy is also supported by 
comparing the red and green curves in Fig. |12[ which cor¬ 
respond to the m = 0 and m = 200 states, respectively. 
The tunneling is well-represented by a function only of 
the transverse mode energy, even when the nature of the 
excitation (azimuthal or radial) is very different. 

The red and green curves in Fig. correspond to 
Konst = 88 Er, giving a transverse confinement har¬ 
monic oscillator frequency of ^ 400 Hz. In this case, 
the tunneling changes by nearly an order of magnitude in 
the given energy range, which can amount to a significant 
thermal dependence of the effective tunneling amplitude. 
As a point of comparison, for strontium in a magic lat¬ 
tice, 10 Er corresponds to a temperature of « 1.66/iK, 
which is comparable to the operating temperatures of op¬ 
tical lattice clocks [nun]. The blue and magenta curves, 
which correspond to the m = 0 and m = 200 states, are 
computed for Konst = 888 Er, giving a transverse con¬ 
finement harmonic oscillator frequency of ~ 850 Hz. The 
dependence of the tunneling on transverse mode energy 
in the same energy range is now significantly smaller. 
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FIG. 13: Interaction matrix elements in a non-separable op¬ 
tical lattice The s- and p-wave interaction integrals Eqs. ( |14[ |- 
(15 I, for m = 0 (top row) and m = 200 (bottom row), n and 
n' are eigenstate indices, which correlate to transverse mode 
quantnm nnmbers. s-wave (p-wave) interactions show a slow 
decay (growth) with transverse mode energy. 



FIG. 14: Non-separability of ID optical lattice eigenstates 
The von Nenmann entropy of non-separability S and the ge¬ 
ometric non-separability f as a function of Wannier state en¬ 
ergy for Konst = SSEr. Red and blue (green and magenta) 
curves correspond to states of azimuthal quantum number 
m = 0 (m = 200). As was the case for the tunneling (Fig. |12| , 
the non-separability is dominantly a function of transverse 
mode energy. 


less than a factor of two. Hence, the axial motion, repre¬ 
sented through the tunneling properties of the lattice, can 
be effectively decoupled from the transverse motion by in¬ 
creasing the effective transverse confinement frequency. 

In Fig. 13 we show the s- and p-wave integrals 
Eqs. (14|-(15). We use units with (. = 0 (\) 

for s- {p-) wave interactions, where Uho is the harmonic 
length associated with the radial direction, in order to 
facilitate comparison with previous works [in EH] using 
a harmonic oscillator approximation. Similar to the case 
using the harmonic approximation, we find a slow de¬ 
cay of s-wave interactions with mode energy, and a slow 
growth of p-wave interactions with increasing mode en¬ 
ergy. When thermally averaged, this weak energy depen¬ 
dence leads to quantitatively similar behavior between 
the Gaussian, non-separable trap case and its harmonic 
oscillator approximation when the temperature is com¬ 
parable or larger than the transverse mode spacing. 


C. Non-separability and anharmonicity of optical 
lattice states 

We now consider the quantitative non-separability of 
the eigenstates of Eq. (§. We decompose a band of states 
as 

V'm,.(r)=5]4'")(r)4-)(z) . (53) 

n 


for the tunneling, the non-separability is dominantly a 
function of the transverse mode energy, irrespective of 
its character (azimuthal or radial), as can be seen by 
comparing the red and blue curves (corresponding to m = 
0) with the green and magenta curves (m = 200). 

To quantitatively assess the degree of anharmonicity 
of the optical lattice eigenstates, we compute the overlap 
of the radial part of these states with harmonic oscilla¬ 
tor functions chosen to match the local curvature of the 
potential. Namely, we compute 

m(i) = 

i/,nho / ^ ^y,n,q,nz 

TLz ,71 


drTZ„{r)(l)ni,, (r) 


(54) 


in the notation of Eq. ( |47| ). In the lowest-order harmonic 
expansion of the potential, ■ The function 

is shown in Fig. for Konst = 88Eij, where 
n is the approximate radial quantum number. The val¬ 
ues at quasimomentum q = —O-Grr/a are representative 
of general quasimomentum away from high-symmetry 
points. For low-lying eigenstates, « ^n,nho: 

showing the near-harmonic nature of these states. How¬ 
ever, in higher-lying eigenstates, the character of the 
wavefunction is spread over many harmonic oscillator 
modes, demonstrating signihcant anharmonicity. 


VI. CONCLUSIONS AND OUTLOOK 


The non-separability measures for this decomposition are 
shown as a function of Wannier state energy for Konst = 
88 Eh in Fig. While the states near the bottom of the 
trap are nearly separable, more highly excited states can 
be very significantly non-separable. As was also found 


We have studied the properties of particles confined in 
non-separable 3D optical potentials, focusing in particu¬ 
lar on a double-well optical tweezer array and a ID op¬ 
tical lattice with transverse Gaussian confinement. Our 
main methodology was discrete variable representations 
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FIG. 15: Anharmonicity of non-separable optical lattice 

states. The overlap of the radial part of the wavefunction, in¬ 
dexed by n, with the approximate harmonic oscillator wave- 
functions, indexed by riho, shows that the radial wavefunc¬ 
tion is spread across many harmonic oscillator modes in ex¬ 
cited states. The values at quasimomentum q = —O.Brr/a are 
representative of general quasimomentum away from high- 
symmetry points. 


separability that quantifies the distance to the nearest 
separable state. Finally, we also presented a variational 
algorithm for determining the nearest separable state to a 
given state, and showed how this nearest separable state 
can be used to gain intuition about anharmonicity, to 
classify quantum states, and to assess the convergence 
of algorithms. We found that the low-lying states of 
the non-separable potentials we considered were nearly 
separable despite often being significantly anharmonic. 
This observation strongly motivates the use of the MPS 
canonical form for non-separable states as a ansatz that 
can be variationally optimized at fixed degree of non¬ 
separability with significantly reduced computational re¬ 
sources compared to direct diagonalization. 

We would like to acknowledge useful discussions with 
Michael Foss-Feig, Adam Kaufman, Brian Lester, Cindy 
Regal, and Xibo Zhang. This work was supported 
by JILA-NSF-PFC-1125844, NSF-PIF-1211914, ARO, 
AFOSR, and AFOSR-MURI. MLW acknowledges sup¬ 
port from the NRC postdoctoral fellowship program. 


(DVRs), which couple the rapid convergence of spectral 
methods with the flexibility and simplicity of grid-based 
methods, as well as variational methods built on top of 
DVR approaches. We found that parameters relevant 
to the construction of effective many-body models, such 
as tunneling and interaction matrix elements, can be 
significantly different from their separable counterparts. 
In particular, we found that lowest-band tunneling am¬ 
plitudes in a non-separable lattice can be increased by 
nearly an order of magnitude for a state with signifi¬ 
cant transverse mode energy compared to a state in the 
transverse ground state. Similarly, we found the lowest- 
lying state with an excitation transverse to the tunnel¬ 
ing direction in a double-well optical tweezer has a tun¬ 
neling rate « 20% larger than the ground state. The 
fact that tunneling depends on transverse motional state 
could have a range of applications, such as thermome¬ 
try and quantum simulation of multi-component systems 
with effective mass imbalance. Interactions were found 
to be less sensitive to non-separability and anharmonicity 
compared with tunneling amplitudes. 

In addition to discussing how effective model parame¬ 
ters change with the trap variables, we also presented a 
quantitative analysis of the non-separability of individ¬ 
ual eigenfunctions by adapting methods from the the¬ 
ory of matrix product states (MPSs). In particular, we 
developed a canonical form for non-separable states in 
terms of a contraction of low-rank tensors describing mo¬ 
tion along each independent direction, and discussed how 
this canonical MPS form is useful for storage, computa¬ 
tion, visualization, and quantification of non-separability. 
Based on this canonical form, we discussed three mea¬ 
sures of non-separability: the Schmidt rank and von 
Neumann entropy of non-separability, both of which de¬ 
pend on a specific bipartition of degrees of freedom, and 
the geometric non-separability, a global measure of non- 


Appendix A: Interaction parameters for radial 
functions 

As mentioned in Sec. |III A~^ the Bessel DVR grids are 
different for different angular momentum m, and so there 
is no single DVR quadrature which can be used to inte¬ 
grate products of functions which have different values of 
m. Hence, in order to find interaction matrix elements 
between radial functions, we must use some other basis 
set. A natural choice is to use the sine DVR in r, adapted 
for odd-parity potentials to ensure that all radial func¬ 
tions vanish at r = 0. Using the sine DVR, all radial 
functions are expressed in terms of the same grid, and so 
the sine DVR quadrature can be used to find interaction 
matrix elements. This procedure works well for \m\ > 3, 
but is slowly convergent for \m\ < 2. This is due to the 
fact that the sine DVR basis functions do not have the 
appropriate asymptotic behavior near r = 0. As |m| in¬ 
creases, the wave function is very small near this region 
due to the centrifugal potential, and so the mismatch of 
boundary conditions does not affect the convergence of 
the algorithm. For the states with \m\ < 2, we use the 
following alternate procedure: 

1. Solve for the eigenstates of the Gaussian potential 
using the Bessel DVR. 

2. Solve for the eigenstates of a radial harmonic oscil¬ 
lator with the same local curvature near the poten¬ 
tial minimum as the Gaussian potential using the 
Bessel DVR. 

3. Find the expansion of the Gaussian potential states 
in terms of the harmonic oscillator states using the 
Bessel DVR quadrature. 

4. Solve for the eigenstates of a Gartesian harmonic 
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oscillator with the same frequency as the radial os¬ 
cillator using the sine DVR. 


5. Use the Cartesian harmonic oscillator functions to 
obtain the radial harmonic oscillator functions on 
an equally spaced radial grid using transformation 
Eq. (lAll) below. 


6. Using the results of step 3 and 5, find the values of 
the Gaussian potential states on an equally spaced 
radial grid. 


One may object that steps 3-5 are unnecessary, as the ex¬ 
act wave functions for the harmonic oscillator are known. 
However, evaluation of these wavefunctions directly in 
terms of orthogonal polynomials is numerically unstable. 
In contrast, using DVR-based eigenvalue methods is nu¬ 
merically stable, even for very highly excited states. 

The expansion of the radial harmonic oscillator states 
\m,nr) in terms of Cartesian harmonic oscillator states 
\p)\q), where “p” labels the x harmonic oscillator state 
and “q” labels the y harmonic oscillator state may be 
accomplished as 


n,, + |m|/2 

\m,nr)= 

/c=—n^ —|m|/2 


m}2,k 



+ k)\nr + 



(Al) 


k) , 


where (9) is the Wigner little-d matrix [37]. Since 
we do not need the full 2D wavefunction but only the 
radial function, we can consider |m, Ur) along the line of 
polar angle (/) = 0 in which r = x. Here, we have 


nr-\-\m\/2 

(r, </, = 0|m, nr) = Y. (-!)"’■ 

fc ———| m |/2 

(A2) 

^ (“0 ^ A:)(0|n^ + ^ - fc) • 


Due to the fact that {x\n) is an odd (even) function if n 
is odd (even), only k such that (nr + ^ — k) is even con¬ 
tribute to the sum, and hence the radial wavefunction is 
real. In contrast to the recurrence relations required for 
the harmonic oscillator wavefunctions themselves, d ma¬ 
trices have a numerically stable recurrence relation [38j . 

At the end of this procedure, we have the values of 
the radial functions evaluated on a grid of equally spaced 
radial points, and so we could think of this as being an 
expansion of the radial wave functions in terms of an odd- 
parity sine DVR basis and use the associated quadrature 
to evaluate overlaps and derivatives. However, due to 
the mismatch in boundary conditions between the ra¬ 
dial function and the sine DVR functions, such an ex¬ 
pansion is inaccurate, and instead standard grid-based 
techniques for integration, e.g. Simpson’s rule, and eval¬ 
uation of derivatives by high-order finite differencing will 


yield more accurate results. In contrast to performing 
derivatives and integration with DVR quadrature, where 
the calculation converges exponentially fast, this proce¬ 
dure features only algebraic convergence with the step 
size Ax, with the particular rate of convergence set by 
the finite differencing and integration scheme. 


Appendix B: Obtaining the MPS form of a 
non-separable state and a variational algorithm for 
finding the nearest separable state 


Similar to the case of the Schmidt form, the MPS rep¬ 
resentation of a quantum state can be obtained via the 
singular value decomposition (SVD). In particular, let us 
assume we have a state ipxi,yj,zt,: where Xi, yj, and Zk 
are discrete indices running over some finite set of basis 
functions (e.g. DVR functions). Then, we can obtain the 


representation Eq. (|M|) as 




(Bl) 

Zfly,, {Zk) 

= '^liyz,Xk ^ = U(xi,yj),vSv 

(B2) 

Axi,{yj,v) 


(B3) 

XkL,cy i^i) 

^Xijflxy 5 ^k'xy^f-l-xy jiVj 



(B4) 


Here, (a, b) denotes the Kronecker product of the indices 
a and b, and SVI^ denotes matrix decomposition of the 
left hand side into the right hand side via the SVD. The 
particular form of Eq. ( |3^ obtained with Eqs. (B1 )-(B4) 
is a mixed canonical form [12] with the gauge conditions 
J dxX^ (^) (^) — J" dzZ^ (z) Zi, (z) — and 

= 1- 

Finding the nearest separable state for a multi-partite 
system is not as simple as for the two-partite case [3(71 l3T] . 
However, given the MPS form of a non-separable state 
Eq. (39), a variational algorithm for obtaining the near¬ 
est separable state (r|'(/'sep) = X {x)Y (y) Z (z) can be 


devised by optimizing the tensors X, Y, and Z individu¬ 
ally in a round-robin fashion. Such an alternating least- 
squares algorithm is similar to the local energy optimiza¬ 
tion coupled with sweeping over lattice sites used in the 
DMRG algorithm of condensed matter physics [T2|. The 
optimal local tensor updates in the case at hand are 


A (a.) = ^ X,^^ (x) ■ y) [z,^^ ■ Z) , 

(B5) 

Y{y)= Y {y) [Z^^^ ■ z) , 

k^xyk'yz 

(B6) 

Z{z)= Y 

l-^x y l-^y z 

(B7) 
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where the updated tensor, e.g., X (x) is normalized af¬ 
ter the update and A • A is shorthand for f d^A(^)A(^). 
Convergence can be assessed by stationarity of the func¬ 
tional llV'sep) — IV')I^- An appropriate starting guess is 
X(x) = Xi(x), Y(y) = Z (z) = 


Z\ (z), which would be the expectation for the nearest 
separable state based on applying the optimal truncation 
at each bipartition. This algorithm can also be straight¬ 
forwardly generalized to find the nearest state with fixed 
Schmidt rank non-separability Xxy and Xyz- 
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